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We present the virial equation of state of low-density nuclear matter composed of neutrons, pro- 
tons and alpha particles. The virial equation of state is model-independent, and therefore sets a 
benchmark for all nuclear equations of state at low densities. We calculate the second virial coef- 
ficients for nucleon-nucleon, nucleon-alpha and alpha-alpha interactions directly from the relevant 
binding energies and scattering phase shifts. The virial approach systematically takes into account 
contributions from bound nuclei and the resonant continuum, and consequently provides a frame- 
work to include strong-interaction corrections to nuclear statistical equilibrium models. The virial 
coefficients are used to make model-independent predictions for a variety of properties of nuclear 
matter over a range of densities, temperatures and compositions. Our results provide constraints on 
the physics of the neutrinosphere in supernovae. The resulting alpha particle concentration differs 
from all equations of state currently used in supernova simulations. Finally, the virial equation of 
state greatly improves our conceptual understanding of low-density nuclear matter. 
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I. INTRODUCTION 

What do we mean by nuclear matter at subsaturation 
density? This is an important conceptual question. The 
binding energy of uniform nuclear matter has a minimum 
at saturation density rig rj 0.16 fm~ 3 . Below this density, 
the system can minimize its energy by forming clusters. 
Therefore the study of low-density nuclear matter is a 
study of clusters. What are the properties of these clus- 
ters? What role do they play in astrophysics? What is 
the relationship between clusters formed in the thermo- 
dynamic limit for an infinite system and those formed in 
laboratory heavy-ion collisions? In this paper, we answer 
these questions using the virial expansion for low-density 
nuclear matter. 

Simulations of core-collapse supernovae 0, Q, giant 
stellar explosions, depend on properties of low-density 
nuclear matter near the neutrinosphere. If one views a 
star in visible light, one sees the photosphere. This is 
the surface of last scattering for photons. Supernovae 
emit 99% of their energy in neutrinos, and if one views 
a supernova in neutrinos, one sees the neutrinosphere. 
This is the surface of last scattering for neutrinos and 
occurs at a density where the neutrino mean-free path 
is comparable to the size of the system. Direct evi- 
dence for the temperature of the neutrinosphere yields 
roughly T w 4 MeV, from the small number of neu- 
trinos detected in SN1987a 0, 0. The neutrinosphere 
density, n ~ 10 11 g/cm 3 ~ 1/1000 no, then follows from 
known cross sections for neutrinos with these energies. 
In this paper, we use the virial expansion to make model- 
independent predictions for the equation of state and 
other properties of low-density nuclear matter near the 
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neutrinosphere. We will use these results in later works 
to calculate neutrino interactions. 

Nuclear statistical equilibrium (NSE) models are com- 
monly used in nuclear astrophysics |5| . These models de- 
scribe low-density nuclear matter as a system of non- or 
minimally-interacting nuclei with a distribution of neu- 
tron number N and charge Z determined from nuclear 
masses through statistical equilibrium. However, NSE 
models fail in an uncontrollable manner as the density 
increases and strong interactions between nuclei become 
important. One would like a model- independent way to 
include strong interactions and a well-defined criterion 
to determine the range of validity of NSE models. In 
this paper we take a step towards this goal by consid- 
ering a simple NSE model with neutrons (n), protons 
(p) and alpha (a) particles. The virial expansion is then 
used to systematically incorporate strong interactions us- 
ing nucleon-nucleon (NN), Na and aa elastic scattering 
phase shifts. 

Pure neutron matter at low density is a very inter- 
esting nearly-universal Fermi system. In the limit that 
the neutron-neutron scattering length is large compared 
to the interparticle spacing and the range of the inter- 
action is small, the properties of dilute Fermi gases are 
expected to be universal, independent of the details of 
the interaction. In a separate paper Q , we have calcu- 
lated the equation of state of low-density neutron matter 
using the virial expansion, in order to assess quantita- 
tively how close the system is to universal behavior at 
finite temperature. The equation of state of resonant 
and dilute Fermi gases has also been studied in labora- 
tory experiments with trapped atoms 0iH,0,^3j an d Ho 
et al. have used the virial expansion to describe Fermi 
gases at hi gh temperatures in the vicinity of Feshbach 
resonances jlllll2|. 

There are many theoretical approaches to low-density 
nuclear matter. Microscopic calculations start from 
NN and three-nucleon interactions that reproduce NN 
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scattering and selected few-nucleon data. In addition 
to the conventional variational [l3l Il4j and Brueck- 
ner 0,0] calculations, renormalization-group methods 
coupled with effective field theory provide new insights 
to nuclear forces and the possibility of a perturbative and 
thus systematic approach to nuclear matter with theoret- 
ical error estimates However, all wave functions 
that are typically used only include two-nucleon corre- 
lations and omit four-nucleon correlations that could be 
important to describe alpha particle formation at low 
densities. We note that the equation of state calculation 
of Buchler and Coon [l^ is close in strategy to the virial 
equation of state, but it takes into account Pauli blocking 
on the phase shifts and thus uses a model NN interaction. 
Moreover, there are a variety of Skyrme-type USUI HI 
and relativistic mean-field |23j , |24{,|25{ parametrizations of 
the nuclear energy functional, which are used to calculate 
ground-state energies and densities of intermediate-mass 
and heavy nuclei. However, it is not clear what form to 
take for the energy functional for a gas of nucleons at 
low densities, where the system forms clusters and the 
binding energy per nucleon is density-independent with 
E/A w -16MeV. 

There are also some phenomenological models that at- 
tempt to provide a reasonable behavior for the equation 
of state over a wide range of densities, temperatures and 
compositions. The Lattimer-Swesty equation of state j2(J 
is based on an extended liquid-drop model. This equa- 
tion of state is almost universally used in modern su- 
pernova simulations Q, Q . In addition, the Shen-Toki- 
Oyamatsu-Sumiyoshi equation of state j2?J is based on 
an approximate Thomas-Fermi calculation using a rela- 
tivistic mean-field interaction. However, it is not clear 
how reliable these models are, and there does not appear 
to be any way to systematically improve them. 

Recently, there are lattice simulations for neutron mat- 
ter using effective field theory j2^|, where the couplings 
are fitted to reproduce NN scattering regularized on the 
same lattice. These calculations are promising, but are 
presently limited to low orders in the effective field theory 
expansion and to small lattices. It is very useful to com- 
pare these lattice results to our virial equation of state 
at low densities. In the virial approach it is simple to 
assess how errors in the NN phase shifts, because of a 
low-order truncation in the effective field theory, impact 
the equation of state. The virial equation of state thus 
provides a valuable check for the lattice results. 

The virial expansion is a very general and old method 
for determining the equation of state of a dilute gas. 
There are two major assumptions. First, that the sys- 
tem is in a gas phase and has undergone no phase tran- 
sition with decreasing temperature or increasing density. 
Second, that the fugacity z = e^l T is small, so the par- 
tition function can be expanded in powers of z. Here fj, 
denotes the chemical potential and T is the temperature. 
The second virial coefficient 62 describes the z 2 term in 
the fugacity expansion of the partition function, while 
the third virial coefficient describes the z 3 term. A gen- 



eral formula relating &2 to the two-body elastic scattering 
phase shifts has been known for some time 1291 . Little is 
known about the third virial coefficient [3fll3l| . The first 
application of the virial expansion to hot matter formed 
in heavy-ion experiments was by Pratt et al. |33| . In ad- 
dition, Venugopalan and Prakash [34|. for example, use 
the virial expansion to study matter formed in relativis- 
tic heavy-ion collisions. In general, however, the latter 
work is interested in much higher temperatures, where 
thermal mesons are important. 

A great advantage of the virial expansion is that it 
correctly includes both bound states and scattering reso- 
nances on an equal footing. In contrast, NSE models only 
include bound state contributions. Since the thermody- 
namics of the system is continuous when a bound state 
moves into the resonant continuum, both weakly bound 
states and low-energy scattering resonances are equally 
important. This was explicitly demonstrated for Fermi 
gases in the crossover region of a Feshbach resonance by 
Ho and Mueller The virial expansion in this paper 
includes contributions from the two-nucleon (deuteron) 
and four-nucleon (alpha particle) bound states and the 
relevant low-energy scattering resonances: the dominant 
ones being the 1 So two-nucleon resonance ( 2 He or 2 n), 
the P3/2 Na resonance and the aa + resonance that is 
the ground state of 8 Be. 

This paper is organized as follows. In Section [HI we 
present our virial formalism. Results are presented in 
Section lTTTl for the virial coefficients and equilibrium prop- 
erties, such as the composition, pressure, entropy, energy 
and symmetry energy. We conclude and give an outlook 
in Section llVl 

II. FORMALISM 

In this section, we introduce the virial expansion for 
a system of neutrons, protons, and alpha particles. We 
first discuss our choice of n, p and a particle degrees of 
freedom. Then we present the virial equation of state 
expanded to second order in the fugacities and calculate 
the entropy and the energy. Finally, we relate the second 
virial coefficients to the relevant scattering phase shifts. 

A. Nuclear Statistical Equilibrium 

There is a close relation between the virial expansion 
and NSE models. The latter include a variety of bound 
states, while the virial expansion includes both bound 
states and scattering resonances. The virial equation of 
state can therefore also be regarded as a systematic ex- 
tension of NSE models to take into account strong inter- 
actions in the resonant continuum. 

The deuteron is automatically included in the virial 
expansion as a bound-state contribution to the second 
virial coefficient 62 ~ e Ed ^ T with experimental deuteron 
binding energy Ed — 2.22 MeV. However, the a particle 
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has a large binding energy, E a = 28.3 MeV, and this will 
only be included in the fourth virial coefficient describ- 
ing four-nucleon interactions. In general, fourth-order 
contributions are expected to be small at low densities 
and high temperatures. However, for lower temperatures 
e E a /T j g ver y l ar ge. This will lead to an unnaturally large 
fourth virial coefficient that will greatly reduce the radius 
of convergence of the nucleonic virial expansion. 

Our solution to this problem is to rearrange the orig- 
inal nucleonic virial expansion and include a particles 
explicitly. This allows us to retain the large e Ea l T terms 
at lower order in the virial expansion. This procedure 
can be extended to explicitly include a variety of addi- 
tional heavy nuclei that may be important at even higher 
densities. We discuss the effects of heavier nuclei and of 
three-nucleon contributions in more detail below. In this 
first paper, for simplicity, we limit ourselves to n, p and a 
particle degrees of freedom. We note that three-nucleon 
bound state effects are expected to be small, because 
their binding energies are small compared to E a , and the 
system will thus favor the formation of alpha particles. 

Therefore, we consider a low-density gas of n, p and 
a particles. In chemical equilibrium, the neutron //„, 
proton n p and alpha [i a chemical potentials satisfy 

+ 2/i„ = fx a , (1) 

so that the neutron z n = e^ n ^ T , proton z p = e^vl T anc l 
alpha z a = e (^+ £ °)/ T fugacities are related by 



2 2 E a /T 



(2) 



As discussed above, we will consider z a to be the same or- 
der as z n or Zp and expand the partition function through 
second order in the fugacities z n , z p and z a . This is equiv- 
alent to the virial expansion with only nucleon densities 
at high temperatures and reduces to the virial expansion 
of a pure a gas at low temperatures (assuming equal num- 
bers of n and p, z n = z p ). Moreover, it is equivalent to 
including only the bound-state contribution in the fourth 
nucleon virial coefficient. 



B. Virial Equation of State 

We expand the grand-canonical partition function Q 
for a system of n, p and a particles confined to a volume 
V as 

Q\%ni Zoli Vi ^) — 1 H~ ZnQn H - ZpQp %aQa 
~\~ ^ZpZnQnp ~l~ ^Zn^aQ an 



H~ ^ZpZaQ ap H~ ) ■ 



(3) 



Here, Qi = Qi(V, T) and Q n , Q p and Q a are the partition 
functions for single particle n, p or a systems 



Qn Qp 2 Qa 1 



(4) 



where A (A Q ) denotes the nucleon (a particle) thermal 
wavelength, 



A = {2n/mT) 1/2 and X a = (2tt/to q T) 1 



/2 



(•5) 



We use m n = m p for the nucleon mass m and m a — 
4m for the alpha particle mass, thus X a = A/2. Be- 
low, we will relate the two-particle partition functions 
Qij = Qij{y,T) with i,j — n,p,a to the second virial 
coefficients. 

Next, we expand log Q to second order in the fugacities 



log Q o ^ 



Qnn 2 Qn 



V A 3 A 3 1 ~ n v 



+ z 
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V 



Qan oQaQ n. 
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We define the second neutron virial coefficient b r , as 



(6) 



A 



— OT , {Qnn r,Qn) ~ r,Tr(Qpp r,Qp)' (^) 



2V 



2V 



where the second approximation is due to neglecting 
the Coulomb interaction between protons and assumes 
charge-independent nuclear interactions. Likewise, the 
second alpha virial coefficient b a is given by 



A 3 I 

~~^~{Qaa ~^Q 



2 ) 

a/ 1 



(8) 



and the virial coefficients describing strong interactions 
between pn and Na particles are 



A 3 I 

bpn = ~^y(Qnp ~ ~^QnQp) ■ 

A 3 I 
^ati ~^~{Qa7i ~^QaQn) 



0) 

^r(Qap - T^QaQp) , (10) 



where again the second approximation neglects the 
Coulomb interaction and strong-interaction charge de- 
pendences. Coulomb effects typically depend on the ra- 
tio of the Coulomb energy to the thermal energy I\ = 
Zfa/mT [52, where Zi is the ion charge, a the fine 
structure constant, and a.; = (3/47rni) 1 / 3 with n.; the 
ion density. For example, for symmetric nuclear mat- 
ter (assumed composed only of nucleons) at a density 
n = 10 12 g/cm 3 and T = 4 MeV one has T p = 0.05, 
and thus we expect Coulomb effects to be small. How- 
ever, this should be checked explicitly in future work and 
Coulomb interactions may be more important when we 
extend the virial approach to densities where heavy nu- 
clei are present. Note that the choice of A 3 in the last 
equation (instead of A 3 ) is our convention. Below, we will 
calculate these second virial coefficients from microscopic 
NN, Na and aa elastic scattering phase shifts. 
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The pressure P is given by the logarithm of the par- 
tition function, which we truncate after second order in 
the fugacities. In terms of the second virial coefficients, 
the pressure can be written as 

jT = y = ~)Ji^ Zn Zp ( Zn z p)^ n 2zpZ n bp n ) 

+ T3~( Zq + Z a b a + 2z a (z n + Z p )b an ) . (11) 

The n, p and a densities follow from derivatives of log Q 
with respect to the chemical potential or the fugacity, 



d logg 

dzi V 



(12) 



V,T 



The resulting n, p and a particle densities are given by 
2 ^ 

n n — ^ (z n -\- 2z n 6 n -\- 1zpZ n bp n -\- Sz a z n b an ) , (13) 
o 

(zp 4" 2z p b n -(- *2zpZ n bp n -(- &z a Zpb an ^ , (14) 



p — A 3 ' 



Ait 



(15) 



where we have used A 3 /A 3 = 8. 

These equations provide a parametric form for the 
virial equation of state that is thermodynamically con- 
sistent. For values of z„ and z p , the alpha fugacity z a 
follows from chemical equilibrium z a = z p z\ e E °>/ T and 
one can then calculate the pressure P and composition, 
n n , n p and n a , from Eqs. (|1 111311 5|) . The total baryon 
density rib is given by 



rib = n n + n p + An a , 



(16) 



and the proton fraction Y p , the number of protons per 
baryon, by 



Y p = {rip + 2n a )/n b . 



(17) 



The dependence of the baryon density and the proton 
fraction on z n and z p can be inverted to yield the virial 
equation of state in terms of 

P = P(z n (n b ,Y p ,T), z p (n b ,Yp,T), z a (z T),T). 

(18) 

Note that these dependences are nonlinear and strongly 
temperature dependent, mainly due to the e Ea l T term 
in the chemical equilibrium condition for z a . In practice, 
for a given proton fraction, we use Eq. I|17|) to deter- 
mine the proton fugacity as a function of the neutron one 
z n (zp, Y p ,T), and generate the virial equation of state in 
tabular form for a range of z n . This maintains the ther- 
modynamic consistency of the virial equation of state. 



for spin-zero particles |29j to include spin and isospin. 
The virial coefficient is related to the partition function 
of the two-particles system Estates e~ E2 ^ T , where the 
sum is over all two-particle states of energy E%. This 
sum can be converted to an integral over relative mo- 
mentum k weighted by the density of states. The differ- 
ence between the density of states of an interacting and a 
free two-particle system can be expressed in terms of the 
derivative of the two-body phase shift dS(k)/dk [29l l35| . 
Finally, if one integrates by parts, the virial coefficients 
can be calculated from an integral over the scattering 
phase shifts summed over all partial waves with two-body 
spin S, isospin T, orbital angular momentum L and total 
angular momentum J allowed by spin statistics. 

The virial coefficient b n describes the interactions in 
pure neutron matter, and we have 



b n (T) 



5/2 



(19) 



Here, — 2~ 5 / 2 is the free Fermi gas contribution due to the 
Pauli principle and (S^ ot (E) is the sum of the T = 1 elastic 
scattering phase shifts at laboratory energy E. This sum 
includes a degeneracy factor (2 J +1), 

COEH ^(2J + 1)J 2S+1L ,(^) 

S,L,J 

= Si So + &ip + 3 <5a Pl + 5 S:ip 2 + 5 (5i Da + . . . 

(20) 

Here and in the following we have neglected the effects of 
the mixing parameters due to the tensor force. We expect 
that their contributions to the second virial coefficients 
describing spin-averaged observables vanish. 

The pn virial coefficient can be decomposed as 



b pn (T) = b nuc (T) - b n (T) . 



(21) 



where 6 nuc is the second virial coefficient for symmetric 
nuclear matter. With b pn and b n one can describe asym- 
metric matter with arbitrary proton fraction. For & nuc , 
the contributions from the deuteron bound state and the 
scattering continuum are given by 

b nuc (T) = ^(e E «/ T -l)-2-^ 

i r°° 

+ ¥jhfj dE ^ E/2TS nuc(E), (22) 



C. Virial Coefficients 

We relate the second virial coefficients to the relevant 
scattering phase shifts. This extends the standard results 



where the term —1 in the deuteron contribution comes 
from the partial integration (the phase shift at zero en- 
ergy being ir times the number of bound states) and 
the factor 3 counts the spin-isospin degeneracy of the 
deuteron. The total phase shift S^ C (E) for nuclear mat- 
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ter also receives contributions from T = states, 

Cc(£) = E (2J+ 1)(2T + 1) W Lj (^) 

= 3(5i So + 3<y 3Sl + 3 5i Pl 

+ 3 <*3 Po + 9 Ss Pl + 15 S 3 p 2 + 15 5i U2 

+ 3& Dl +5 5s D2 +7*8 D , + ... (23) 

For spin-zero alpha particles we can directly use the 
results of |2^, [3^ , which give for the virial coefficient b a 

nl/2 roo 

b a (T) = — / dE e- E ' 2T 5^ (E) + 2- 5 /\ (24) 

where +2~ 5 / 2 describes the contribution for a free Bose 
gas and the total phase shift 5 t ° t {E) for elastic aa scat- 
tering is given by 

5 t ° t (E)=J2(2L + l)5 L (E) 

L 

= (53 + 550 + 950 + 13(51 + ... (25) 

For the aa virial coefficient and for the following Na 
virial coefficient, we have not taken into account the ef- 
fects of inelasticities. For low temperatures, we expect 
their contributions to be small due to the tight binding 
of the alpha particle. 

Finally, the Na virial coefficient is given by 

b an (T)= (Jj^^fdEe-^C^E), (26) 

with total phase shift (E) for Na scattering at nucleon 
laboratory energy E, 

L,J 

= 2 s s 1/2 + 2 (5 Pl/2 + 4 Sp 3/2 + 4 <5d 3/2 + 6 <5 D5/2 
+ 6<5 F5/2 +8<5 F7/2 +... (27) 

Once the four virial coefficients b n (T), b pn (T), b a (T) and 
b an (T) have been calculated, the pressure is determined 
from Eq. (|ll|l using fugacities z n and z p that reproduce 
the desired baryon density rib and proton fraction Y p us- 
ing Eas. (12113111 . 

D. Entropy and Energy 

The entropy S and the energy E can be calculated from 
the virial equation of state using thermodynamics [TH 
l35| . The entropy density s = S/V follows from 




Here the temperature derivative is at constant fi p and 
/i„, and thus constant \x a due to chemical equilibrium. 
This leads to 

5P 

s = — - n n log z n - n p log z p - n a log z a 
T 

+ -tt {z 2 a b' a + 2z a (z n + z p )b' an ) , (29) 

"a 

where b\ denotes the temperature derivative of the second 
virial coefficients b'JT) = dbi{T) / dT . The energy density 
e = E/V can be calculated from the entropy density by 

e = Ts + 2J n i^i ~ P 

3 2T 2 
= — P — n a E a H — g- ((z n + z p )b n + 2zpZ„6p n ) 




We emphasize that the a particle binding-energy contri- 
bution n a E a is very important for the energy. Finally, 
the entropy per baryon S/A, energy per baryon E/A and 
the free energy per baryon F/A are given by 

S s E e F f 

-j = — , = — and — = — , (31) 

A rib A ri}, A rib 

with the free energy density / = e — Ts. 

III. RESULTS 

In this section, we present results for the second virial 
coefficients, the equation of state and the composition, 
as well as the entropy, energy and symmetry energy. 

A. Virial Coefficients 

We first calculate the NN virial coefficients. We take 
the NN p hase shifts from the Nijmegen partial wave anal- 
ysis |36( and use pn phase shifts for both T = and 
T = 1 states. We thus neglect the Coulomb interaction 
and strong-interaction charge dependences. In general, 
we expect that their effects on interaction energies are 
small. All partial waves with L ^ 2 are included and we 
have checked that higher contributions are negligible for 
the temperatures of interest. The resulting total phase 
shifts for neutron and nuclear matter are shown in Fig. ^ 

For neutron matter, we observe that (5^ ot (E) is approx- 
imately energy-independent over a wide range. The de- 
crease of the phase shift with increasing energy is 
compensated by the contributions from higher angular 
momenta. As a result, the neutron virial coefficient b n 
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TABLE I: NN virial coefficients for different temperatures. In addition to the full results, we also give the virial coefficients 
calculated only from the large S-wave scattering lengths (ai So = —23.768 fm, <23 Sl = 5.420 fm) and the deuteron binding energy. 
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will be approximately temperature-independent. In con- 
trast, 5^ C (E) decreases with increasing energy due to the 
decrease in the additional 3 Si phase shift. This behav- 
ior of the total phase, along with the strong temperature 
dependence of the deuteron contribution e Ed ^ T , will lead 
to second virial coefficients b pn and 6 nuc that decrease 
rapidly with temperature. 

In Table [IJ we give our results for the NN virial coeffi- 
cients b n (T), b pn (T) and b nuc (T) 7 as well as their deriva- 
tives Tb'^T). As expected, the NN virial coefficients are 
dominated by the large S-wave scattering lengths and 
deuteron physics, but effective range and higher partial 
wave contributions are significant even for these relatively 
low temperatures. We also emphasize that the effects of 
the low-energy 1 So resonance encoded in b n become more 
important in neutron-rich matter. 

Next, we calculate the Na virial coefficient. We ne- 
glect Coulomb interactions and thus use neutron-alpha 
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FIG. 1: (Color online) The total phase shift <5* ot (£) for neu- 
tron matter and 5^ C (E) for nuclear matter versus laboratory 
energy E. For reference, we also show the contributions from 
only the S-wave phase shifts. 



phase shifts up to L ^ 3. For laboratory energies 
E < 20 MeV, we take the simple phase shift fits of Arndt 
and Roper [3jj , which are constrained by scattering data 
over these energies. We extrapolate to higher energies us- 
ing microscopic phase shift predictions based on optical 
model calculations by Amos and Karataglidis .38]. Note 
for these calculations we assume |<5l. 7 | < 7T. Since we 
neglect the effects of inelasticities, we take the real part 
of the predicted phase shifts of Amos and Karataglidis. 
The differences between the real part and the absolute 
value of the phase shifts are generally small over the 
relevant energies. The neutron-alpha elastic scattering 
phase shifts are shown in Fig. [5J The P 3 / 2 wave has 
a resonance at E w lMeV with P3/2 scattering length 
ap 3/2 = — 62.95 fm 3 39], and therefore will be the most 
important contribution to the Na virial coefficient. The 
other partial waves are not resonant at low energies. Our 
results for the Na virial coefficients are listed in Table ITU 
We expect that more accurate high-energy phases will 
only slightly change our results for b an for temperatures 
T > 5 MeV. 

Finally, we show the phase shifts for elastic act scatter- 
ing in Fig. The low-energy phase shifts are taken from 
Afzal et al. j4fj and the phase shifts for energies between 
30 MeV and 70 MeV are from Bacher et al. 41]. The 
phase shifts display a pronounced + resonance at very 
low energies that corresponds to the 8 Be ground state. 
This resonance is crucial for He burning in red giant 
stars. Just below 35 MeV there are also two very close 2 + 
resonances (they are both at E w 34 MeV in Fig. |3Jl and 
in general the phase shifts become large at high energies. 
We include all aa phases with L ^ 6 for the calcula- 
tion of b a , where we make the following approximations. 
First, we neglect the Coulomb interaction. We also ne- 
glect inelastic channels. Fortunately, the a particle is so 
tightly bound that important inelastic contributions do 
not arise until relatively high energies. Finally, we trun- 
cate the integration for b a at 70 MeV (the extent of the 
data) . Our results for the aa virial coefficients are given 
in Table ITTI As for the NN virial coefficients, b an and b a 
are dominated by the relevant low-energy resonances, but 
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TABLE II: Nq and aa virial coefficients for different temperatures. In addition to the full results, we also give the virial 
coefficients calculated only from the dominant low-energy resonances (P 3/ / 2 -wave for No, and S- and D-waves for aa). 



T[MeV] 


b an full 


b a „ only P 3/2 -wave 


Tb'an full 


b a full 


b a only S- and D-waves 


Tb'a full 


1 


1.51 


1.73 


1.21 


2.55 


2.55 


1.59 


2 


2.26 


2.48 


0.90 


4.12 


4.00 


2.95 


3 


2.57 


2.78 


0.63 


5.64 


5.07 


4.81 


4 


2.73 


2.91 


0.44 


7.44 


6.02 


7.90 


5 


2.81 


2.97 


0.32 


9.57 


7.01 


11.3 


6 


2.86 


2.99 


0.23 


11.9 


8.05 


14.3 


7 


2.89 


2.99 


0.18 


14.3 


9.08 


16.3 


8 


2.92 


2.98 


0.15 


16.5 


10.0 


17.3 


9 


2.93 


2.96 


0.14 


18.6 


10.9 


17.5 


10 


2.95 


2.93 


0.13 


20.4 


11.7 


17.0 




E [MeV] E [ MeV l 



FIG. 2: (Color online) The phase shifts for elastic neutron- 
alpha scattering Slj(E) versus laboratory energy E. As 
discussed in the text, the solid lines are from Arndt 
and Roper [23 and the symbols are from Amos and 
Karataglidis |38j| . For clarity, we do not show the F-waves 
included in our results for ban- 



FIG. 3: (Color online) The phase shifts for elastic aa scat- 
tering Sl(E) versus laboratory energy E. As discussed in the 
text, the phase shifts are taken from Afzal et al. [4(1 and from 
Bacher et al. |41| . 



b a receives important contributions from higher angular 
momentum resonances for T > 5 MeV. 

In addition to the values for the NN, Na and aa virial 
coefficients given in Tables HI and ITT1 for T ^ 10 MeV, we 
extend our results to higher temperatures in Table II I II 
We list these results in a separate table, as our approxi- 
mations, mainly the neglect of inelastic channels and the 
truncation of the integration over the phase shifts at the 
extent of the data, may be more severe. 

B. Composition 

We now discuss the a particle concentration. The a 
mass fraction is given by x a — An a /ni,. In Fig.^J we show 



x a for a temperature T = 4 MeV and proton fraction 
Y p = 1/2 with various values of the virial coefficients. 
The largest a particle concentration results for a free n, p, 
a gas, where all bi = 0. Including nucleonic interactions 
b n ,b P n,b nuc 5^ 0, but still keeping b a ,b an — reduces 
x a significantly. This is because the attractive nuclear 
interactions reduce the nucleon chemical potential, and 
this leads to a reduction of x a . Including also b an ^ 
increases the a mass fraction, but x a is still smaller 
than its free value. The latter effect is dominated by 
the low-energy Na resonance. Finally, including b a ^ 
only increases x a by a very small amount. Therefore, in 
general, NN virial coefficients are expected to be most 
important, while the Na virial coefficient leads to small 
changes, and b a is least important. This hierarchy is 
because typically 

The LS model excludes nucleons from the volume occu- 
pied by alpha particles. Although this is a prescription, 
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TABLE III: NN, Na and aa virial coefficients for higher temperatures. Note that inelasticities are neglected and that the 
integration over the aa elastic scattering phase shifts is truncated at E ^ 70MeV. 



T[MeV] 


b n full 


Tb' n full 


b pn full 


Tb' pn full 


6„ uc full 


Tb' nuc full 


b an full 


Tb' an full 


b a full 


Tb' a full 


12 


0.313 


0.013 


1.76 


-0.73 


2.08 


-0.72 


2.97 


0.12 


23.3 


14.7 


14 


0.315 


0.014 


1.66 


-0.65 


1.97 


-0.64 


2.98 


0.10 


25.4 


11.7 


16 


0.317 


0.014 


1.57 


-0.59 


1.89 


-0.58 


3.00 


0.07 


26.7 


8.48 


18 


0.319 


0.013 


1.51 


-0.55 


1.82 


-0.54 


3.00 


0.02 


27.5 


5.44 


20 


0.320 


0.011 


1.45 


-0.52 


1.77 


-0.51 


3.00 


-0.04 


28.0 


2.69 



0-5 1 i- 




n b [fm -3 ] 

FIG. 4: (Color online) The a mass fraction x a versus baryon 
density m, for T — 4MeV, Y p = 1/2 and various values of the 
virial coefficients, as explained in the text. The fugacities are 
small over this density range z„ — z v < 0.05 and z a < 0.004. 



it leads to an a mass fraction that vanishes in the high- 
density limit ^3 • The effects of the excluded volume are 
small at low densities. However, the excluded volume cor- 
responds to a repulsive Na interaction that reduces the 
a particle concentration. In contrast, the neutron- alpha 
resonance and the attractive total phase shift leads to 
an attractive Na infraction in the virial expansion. As a 
result, the a mass fraction increases, as shown in Fig. 0J 
Therefore, the excluded volume alone gives the wrong 
sign for Na interactions in the low-density limit, but we 
caution that the effects of the excluded volume and the 
Na virial coefficient may be small at very low densities. 

(3) 

The third virial coefficient 6„ uc can be used to make 
a simple error estimate of neglected terms in the virial 

(3) 

expansion. From the above conclusion, we expect 6nu C 
to be the most important of the third virial coefficients, 

(3) 

and we therefore take into account only banc i n this es- 
timate. For simplicity, we consider symmetric nuclear 
matter with z p = z n = z, for which the virial equation of 



state, Eqs. (|11I13I15|) . including 6„u C reads 
y = ^(z + z 2 b DUC + z 3 b£> c ) 

+ j^(z a + zlb a + 4z a z b an ) , (32) 
n = n p + n n = ^(z + 2z 2 b nuc + 3z 3 b^ c + 8z a z b an ) , 



n a — — (z a + 2zlb a + iz a zban) ■ (34) 
K 

In Fig. we again plot the a mass fraction. The solid 

(3) 

line is our previous result with &n UC = 0. We also give ap- 
proximate error bands by considering 6 rmc = 10 (which 

gives a lower x a ) and bnuc = —10. This estimates the 
error only from the effect of a typical neglected term in 

the virial expansion. Our estimate |6nuc| < 10 is some- 
what arbitrary, as it is based on the observation that 
all second virial coefficients are comfortably less then 10 
for T = 4 MeV. We see from Fig. that the result- 
ing error band is small. Therefore, our virial expansion 
makes a model-independent prediction for the a particle 
concentration. We also contrast our results in Fig. [S] to 
the x a predicted by two phenomenological equations of 
state. The Lattimer-Swesty (LS) equation of state j2|J 
is based on an extended liquid drop model and is al- 
most universally used in modern supernova simulations. 
The LS model predicts too few alpha particles over the 
densities where the virial equation of state is applica- 
ble. Alternatively, the Shen-Toki-Oyamatsu-Sumiyoshi 
(Shen) equation of state [2?| is based on an approximate 
Thomas Fermi calculation using a relativistic mean-held 
interaction. The Shen equation of state predicts slightly 
too high values for x a at this temperature. 

It is a simple matter to apply the virial equation 
of state to different proton fractions and temperatures. 
In Fig. |SJ we show how the a particle concentration 
decreases with decreasing proton fraction Y p for T = 
4 MeV. In addition, we give x a for Y p = 1/2 and a range 
of temperatures T = 2, 4 and 8 MeV in Fig. □ At the 
lowest temperature T — 2 MeV, x a rises rapidly at low 
densities where the virial coefficients make only small 
contributions. Therefore the error band is very small. 
The error band becomes more important for T = 8 MeV. 
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FIG. 5: (Color online) The a mass fraction x a versus baryon 
density n b for T = 4MeV and Y p = 1/2. The dotted error 
band for the virial equation of state is based on an estimate 
of a neglected third virial coefficient &iu' c = ±10. Also shown 
are results for the LS [2rj| and Shen |27| equations of state. 
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x 




n fa [fm" 3 ] 

FIG. 6: (Color online) The a mass fraction x a versus baryon 
density n b for T = 4MeV and proton fractions Y p = 0.1, 0.2, 
0.3, 0.4 and 0.5. For the lowest proton fraction Y p = 0.1, the 
fugacities are z n < 0.2, z p < 0.005 and z a < 0.001 over this 
density range. 

The LS equation of state is seen to significantly under- 
estimate x a at all three temperatures. This may be due 
to a simple oversight in the published version |2(| of the 
LS equation of state. Lattimer and Swesty include the 
neutron-proton mass difference in the proton chemical 
potential, but do not appear to include twice this dif- 
ference in the alpha chemical potential. In contrast, the 
Shen equation of state is close to our virial results at 
low densities. However, both LS and Shen equations of 
state predict too few a particles at high temperatures 
T > 10 MeV. This can be clearly seen from the a mass 
fraction shown in Fig. H for T = 20 MeV. The differ- 
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FIG. 7: (Color online) The a mass fraction x a versus baryon 
density n b for T = 2, 4 and 8 MeV and Y p = 1/2. The dotted 
error bands for the virial equation of state are based on an 

(3) 

estimate of a neglected third virial coefficient 6n UC = ±10. 
Note that for T = 2 MeV the error band is very small. Also 
shown are results for the LS |2G| and Shen |27| equations of 
state. Over this density range, the fugacities are z„ = z p < 
0.03, 0.10, 0.16 and z a < 0.27, 0.09, 0.03 for T = 2, 4, 8 MeV. 



a 

x 




n b [fm" 3 ] 

FIG. 8: (Color online) The a mass fraction x a versus baryon 
density n b for T = 20 MeV and Y p = 1/2. The dotted error 
bands for the virial equation of state are based on an estimate 
of a neglected third virial coefficient b ^lc = ±10. Also shown 
are results for the LS [26f and Shen |27fl equations of state. 
Over this density range, the fugacities < 0.21 and 

z a < 0.07. 

ences of the LS and Shen equations of state at the higher 
densities shown for T = 20 MeV in Fig. [S] are due to the 
differences in the symmetry energy |43j |. 

Both phenomenological equations of state shown in 
Fig. [7] give a particle concentrations that first increase 
with density and then rapidly decrease. This decrease 
is due to the formation of heavy nuclei. Heavy nuclei 
are not included in our first studies. We will extend the 
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FIG. 9: (Color online) The threshold baryon density rib above 
which the Shen equation of state |27| predicts a heavy nuclei 
mass fraction larger than 10%. This figure is for Y p — 1/2. 



virial formalism to include heavy nuclei in future work. 
Therefore, our virial results are currently limited to lower 
densities, where one does not expect large contributions 
from heavy nuclei. In Fig. 03 we plot the density above 
which the Shen equation of state has a heavy nuclei mass 
fraction larger than 10%. This provides an estimate for 
the range of validity of the virial n, p, a equation of state. 
For this estimate we have considered the Shen instead of 
the LS model due to the LS error in x a . 



C. Pressure 

Next, we present our results for the pressure. In 
Fig. HOI we plot the pressure of symmetric nuclear mat- 
ter for T = 10 MeV. The formation of a particles re- 
duces the total number of particles in the system, and 
this significantly lowers the pressure. For comparison, 
we also show in Fig. 1101 the microscopic Fermi hyper- 
netted-chain (FHNC) equation of state of Friedman and 
Pandharipande 13|. The latter calculation is based on a 
variational wave function with only two-nucleon correla- 
tions. However, the model-independent virial equation of 
state shows that there are large a concentrations in this 
density regime. For example, at the threshold density of 
n b = 0.013 fm~ 3 (from Fig. |§J) we have x a = 0.40. Micro- 
scopic calculations may require wave functions with four- 
nucleon correlations to accurately describe these a par- 
ticle contributions. This suggests that all present micro- 
scopic variational or Brueckner equations of state are in- 
complete in their description of low-density nuclear mat- 
ter. In Fig. 1101 we also give the pressure for densities be- 
yond the threshold density of Fig. [5] in order to compare 
virial and FHNC results, which both omit the formation 
of heavy nuclei. 

The error band in Fig. 1101 is again estimated from a 




n b [fm~ ] 

FIG. 10: (Color online) The pressure P of symmetric nuclear 
matter versus baryon density nb for T = 10 MeV. The dotted 
error band for the virial equation of state is based on an es- 
timate of a neglected third virial coefficient b£u C = ±10. Also 
shown are virial results without a particles (z a — 0) and 
the microscopic FHNC results of Friedman and Pandhari- 
pande (FP) 13]. The effects of interactions are illustrated 
by comparing to the pressure of a free gas of nucleons and 
alpha particles. Over this density range, the fugacities are 
z n = z p < 0.21 (< 0.51 without alphas, < 0.32 for the free 
gas) and z a < 0.04 (< 0.17 free gas). 



neglected third virial coefficient &n UC = ±10. For the 

(3) 

pressure, this involves a cancellation. Increasing 6nuc de- 
creases the nucleonic contribution to the pressure due to 
the additional attractive interaction between nucleons. 
However, a larger 6nuc also reduces the formation of a 
particles (due to a lower alpha chemical potential), and 

this acts to increase the total pressure. Therefore, a pos- 
ts) 

itive 6nuc reduces the pressure at very low densities, but 
then leads to a larger pressure with increasing density. 

In Fig. 1111 we show our results for the pressure of sym- 
metric nuclear matter for T = 4 MeV in comparison to 
the phenomcnological LS and Shen equations of state. 
The pressure given in Fig. ^] also includes the contri- 
bution of an electron and photon gas as in Appendix C 
of 26] . We find the pressures given by the virial equation 
of state and the phenomenological models agree well. 

In a separate paper |6J] , we have demonstrated that the 
equation of state of low-density neutron matter scales to 
a very good approximation. Table [I] shows that the neu- 
tron virial coefficient b n (T) is practically independent of 
the temperature. If all virial coefficients are tempera- 
ture independent, and provided the condition for chemi- 
cal equilibrium does not introduce a strong temperature 
dependence through z a (z n , z p , T), then the power series 
in the fugacity for the pressure and the baryon density 
will have no explicit temperature dependence. In this 
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FIG. 11: (Color online) The pressure P of symmetric nuclear 
matter versus baryon density nj, for T — 4 MeV. As discussed 
in the text, these results include the contribution of an elec- 
tron and photon gas, and we also show the virial pressure 
without this contribution. In addition, we give the results of 
the LS |26j and Shen [2jJ equations of state. Over this density 



range, the fugacities are z n 



0.02 



0.015 



< 0.06 and z a < 0.01. 



> 0.01 



0.005 - 




0.005 
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FIG. 12: (Color online) The pressure P of pure neutron mat- 
ter and of symmetric nuclear matter calculated for T = 1, 4 
and 8 MeV and scaled to a temperature of J" = 4 MeV (for 
details see text). The circle indicates where z n = 0.5 for pure 
neutron matter. For symmetric nuclear matter, the fugacities 
are z n = z v < 0.10 and z a < 0.09 over this density range. 



FIG. 13: (Color online) The entropy per baryon S/A versus 
baryon density rib for T = 4 MeV and Y p = 1/2. As discussed 
in the text, these results include the contribution of an elec- 
tron and photon gas. Also shown are results for the LS |2?i |. 
Shen [221 and the Wolff equations of state. The effects 
of interactions are illustrated by comparing to the entropy of 
a free gas of nucleons and alpha particles. Over this density 
range, the fugacities are z n = z v < 0.11 (< 0.17 for the free 
gas) and z a < 0.17 (< 0.81 free gas). 



P(n b ,T')= (TLY P(n b ,T), 



5/2 



(36) 



with n b = (T'/Tf/ 2 n b . 

We study this scaling symmetry in Fig. ^1 where we 
have used Eq. (|36|l to predict the pressure of neutron 
matter for T" = 4 MeV from virial pressures calculated 
for T = 1 and 8 MeV. The agreement with the unsealed 
T = 4 MeV virial pressure is excellent. This demon- 
strates that low-density neutron matter scales, even when 
interactions are included. Fig. ^] also shows scaled pres- 
sures of symmetric nuclear matter for T' = 4 MeV using 
T = 1 and 8 MeV results as input. These scaled pres- 
sures do not agree with the T — 4 MeV virial pressure. 
Therefore, nuclear matter does not scale, which is due 
to clustering and mainly comes as a result of the strong 
temperature dependence of z a and x a . 



D. Entropy and Energy 



case, the pressure will scale as 

PK,T) = T 5 / 2 /K/T 3 / 2 ), (35) 

where f(x) is some function of n b /T z / 2 . If this scaling 
relation holds, one can predict the pressure P(n b ,T') at 
a new temperature T 1 from the original P(n b , T) through 



The entropy per baryon S/A of symmetric nuclear mat- 
ter is shown for T = 4 MeV in Fig. El Here, our virial 
calculation again includes the entropy of an electron and 
photon gas as in Appendix C of |26(. The entropy is 
clearly seen to reflect the composition. At low densities, 
the LS entropy is above the other calculations because 
of the low a mass fraction in the LS model. At high 
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FIG. 14: (Color online) The entropy per baryon S/A versus 
baryon density rib for T = 4 MeV and Y p = 0.05. As discussed 
in the text, these results include the contribution of an elec- 
tron and phot on g as. Also shown are results for the LS |2fjj | 
and the Shen [27| equations of state. The circle indicates 
where z n = 0.5 (z p = 0.002 and z a = 0.001). 
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FIG. 15: (Color online) The energy per baryon E/A versus 
baryon density n b for T = 4 MeV (and also T = 2 MeV (dot- 
ted) for nuclear matter). In order to clearly show the interac- 
tion effects we have subtracted the free kinetic energy 3T/2 
from all curves. The circle indicates where z„ = 0.5 for pure 
neutron matter. For symmetric nuclear matter, the fugacities 
are z„ = z p < 0.08, 0.02 and z a < 0.03, 0.11 for T = 4 and 
2 MeV respectively. 



densities, the entropy in the phenomenological models is 
below our virial result due to the formation of heavy nu- 
clei. This deviation coincides with the threshold density 
of Fig- EI ln Fig- El we also show the entropy per baryon 
for the Wolff equation of state 44], to give a measure 
for the model dependence of the onset of heavy nuclei. 
The entropy per baryon for extremely neutron-rich mat- 
ter Y p — 0.05 is shown in Fig. El Now there is very 
little change in composition, and thus good agreement 
between the virial and phenomenological results. 

In addition to the entropy, we show our results for the 
energy per baryon E/A in Fig. El In the limit of very 
low density, the energy reduces to the kinetic energy of an 
ideal classical gas 3T/2. Therefore, we have subtracted 
3T/2 from E/A in Fig. El i n order to clearly show the 
interaction effects. We find that the energy per baryon 
for pure neutron matter is somewhat below 3T/2 due to 
the attractive nuclear interactions. Our neutron matter 
results are in excellent agreement with the scaling law for 
the energy density e = 3/2P, which follows from Eq. I|30|l 
when the virial coefficients are temperature independent, 
as is the case for neutron matter with b' n « 0. In contrast, 
E/A for symmetric nuclear matter drops very rapidly 
with increasing density. This reflects the large binding 
energy as a particles form. As can be seen from Fig. 1151 
this drop becomes more rapid for lower temperatures. 
We therefore expect that the inclusion of heavy nuclei, 
or many-nucleon correlations, are necessary to obtain a 
constant energy per particle E/A w —16 MeV for T = 
nuclear matter at subsaturation densities. 



E. Symmetry Energy 

The symmetry energy Se characterizes how the energy 
rises as one moves away from equal numbers of neutrons 
and protons, 



1 d 2 E 



dY 2 A 



(37) 



Y p = l/2 



We evaluate the second derivative of the energy per 
baryon with respect to the proton fraction numerically. 
Fig. El shows the symmetry energy for T — 4 MeV. At 
very low density, Se rises slowly with density. To low- 
est order in the density and neglecting a particles, the 
symmetry energy is given by 

S E « | TX 3 n b {b pn - 2T&;„/3 - b n + 2T6'„/3) . (38) 

Then, as a particles form, Se rises much faster with den- 
sity. As a result of clustering, the symmetry energy is 
large even at a very small fraction of saturation density. 
This is analogous to the observed behavior for the energy. 

We also show in Fig. 1161 the symmetry free energy Sf 
defined as 



S F 



1 d 2 F 
SdY^A 



(39) 



Y p = l/2 



To lowest order in the density and neglecting a particles, 
the symmetry free energy is given by 



Sf 



T 
1 



1 



TX 3 n b {b p 



(40) 
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FIG. 16: (Color online) The symmetry energy Se and the 
symmetry free energy Sf versus baryon density for T — 
4 MeV. Also shown are virial results without a particles (z a — 
0). Over this density range, the fugacities < 0.09 

(< 0.37 without alphas) and z a < 0.06. 



where the T/2 term is from the entropy of mixing in the 
free energy density. Again, as a particles form, we find 
in Fig. 1161 that Sf increases rapidly with density. 

Finally, the temperature dependence of the symmetry 
energy is shown in Fig. El for T = 2, 4 and 8 MeV. The 
lower the temperature, the more rapidly Se rises with 
density. We conclude that, in the thermodynamic limit, 
the symmetry energy is large at low densities due to clus- 
ter formation. 

Isospin observables in asymmetric heavy-ion collisions, 
such as the N/Z ration of emitted fragments, are often 
analyzed with semiclassical simulations, where nucleons 
move in isospin-dependent mean fields and undergo two- 
nucleon collisions. The isospin-dependent mean-field im- 
plies a certain symmetry energy. The density depen- 
dence of this symmetry energy can be adjusted to im- 
prove agreement with the data, e.g., see [45j . This phe- 
nomenological symmetry energy at low densities is in gen- 
eral much smaller than our virial Se- This is because the 
symmetry energy in the virial expansion includes the con- 
tributions of clusters with high central densities even at 
low baryon density. The relation between Se and the 
symmetry energy of mean-field models for single nuclei 
should therefore be investigated in future work. 

IV. CONCLUSIONS 

We have presented the equation of state of low-density 
nuclear matter using the model-independent virial expan- 
sion. The virial equation of state is thermodynamically 
consistent and was calculated directly from binding en- 
ergies and scattering phase shifts, without the need of 
a nuclear interaction model or nuclear wave functions. 




n b [fm" 3 ] 

FIG. 17: (Color online) The symmetry energy Se versus 
baryon density for T = 2, 4 and 8 MeV. Over this den- 
sity range, the fugacities are z n = z v < 0.02, 0.09, 0.13 and 
z a < 0.18, 0.06, 0.01 for T = 2, 4, 8 MeV. 



Therefore, our results provide a benchmark for all mi- 
croscopic or phenomenological equations of state at low 
densities. The virial equation of state can also be used 
to constrain models that extend to higher densities. In 
some cases, our virial results may be used to constrain 
model parameters or functional forms for model interac- 
tions. In other cases, more extensive modifications of the 
phenomenological equations of state may be necessary in 
order to reproduce the virial low-density limit. 

Tightly bound nuclei can lead to unnaturally large 
bound state contributions to the higher-order virial coef- 
ficients. Therefore, we have considered a gas of nucleons 
and alpha particles to explicitly take into account the 
large a particle binding energy. This leads to a larger 
radius of convergence compared to a nucleonic virial ex- 
pansion. We have calculated the second virial coefficients 
for NN, Na and aa interactions directly from the rele- 
vant binding energies and scattering phase shifts. For the 
temperatures of interest, the virial coefficients are dom- 
inated by the low-energy bound states and resonances, 
but higher partial wave contributions are significant. 

We have found that the second virial coefficient 
for neutron matter b n is approximately temperature- 
independent, which leads to an approximate scaling sym- 
metry of pure neutron matter. The second virial coef- 
ficients for nuclear matter b pn and 6 nuc decrease with 
temperature, and as a result nuclear matter does not 
scale. The second virial coefficient for Na interactions 
b an is well approximated by the resonant P3/2 contribu- 
tion. This holds over a wide range of energies due to 
a cancellation among the additional partial waves. The 
second virial coefficient for aa interactions b a increases 
with temperature, as the aa phase shifts become large at 
higher energies. The inclusion of the Coulomb interac- 
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tion for the virial coefficients, and the proper treatment 
of inelasticities and mixing parameters (due to the tensor 
force) are important topics for future improvements. 

We have used the virial coefficients to make model- 
independent predictions for a variety of properties of nu- 
clear matter over a range of densities, temperatures and 
composition. The resulting a particle concentration dis- 
agrees with all equations of state currently used in su- 
pernova simulations. The contributions from low-energy 
resonances show most prominently in the composition of 
low-density nuclear matter. 

The effects of the second virial coefficients follows a 
natural hierarchy, with the NN virial coefficients being 
more important than the Ncv virial coefficient, and the 
aa virial coefficient leading to very small changes. We 
have used this hierarchy to make simple error estimates 
by studying the effects of a neglected third virial coeffi- 
cient Ofjuc — ±10. The resulting error bands are small. 
Detailed investigations of the effects of omitted higher 
virial coefficients, as well as the effects of heavy nuclei 
will be the topic of separate studies. For a better error 
estimate, it is important to also have a reliable calcu- 
lation of the third virial coefficients, which could, e.g., 
come from a calculation in the effective field theory for 
halo nuclei |4l|. Furthermore, these error estimates can 
also help clarify the domain of validity of the virial ex- 
pansion. 

The physics of nuclear matter is very different from 
neutron matter due to clustering. As the density in- 
creases, a particles form and this leads to a significant 
reduction of the pressure of low-density nuclear matter. 
Similarly, clustering increases the binding energy and it 



reduces the entropy, which reflects the number of par- 
ticles in the system. As the density increases further, 
heavier nuclei and larger clusters form, and a particles 
become less important. In fact, the breakdown of the 
virial expansion for nuclear matter is due to the forma- 
tion of heavy nuclei, which is reached before the nucleon 
fugacities become large. In nuclear matter at T — 0, we 
expect that one has to explicitly include heavy nuclei, or 
many-nucleon correlations, to obtain a constant energy 
per baryon E/A i=s — 16MeV and a constant symmetry 
energy at low densities. These virial results greatly im- 
prove our conceptual understanding of nuclear matter at 
subsaturation densities. 

The virial equation of state can be regarded as a sys- 
tematic extension of NSE models for strong interactions 
in the resonant continuum. Such a comparison will also 
provide valuable insights how to systematically organize 
the virial expansion to include heavy nuclei. In addition, 
topics of future studies include model-independent pre- 
dictions for the neutrino response of low-density nuclear 
matter and a detailed comparison with nuclear lattice 
calculations [47) . 
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